Self-trapping mechanisms in the dynamics of three coupled Bose-Einstein condensates 
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We formulate the dynamics of three coupled Bose-Einstein condensates within a semiclassical 
scenario based on the standard boson coherent states. We compare such a picture with that of Ref. 
and show how our approach entails a simple formulation of the dimeric regime therein studied. 
This allows to recognize the parameters that govern the bifurcation mechanism causing self-trapping, 
and paves the way to the construction of analytic solutions. We present the results of a numerical 
simulation showing how the three-well dynamics has, in general, a cahotic behavior. 
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An increasing interest for the dynamics of coupled 
bosonic wells [known in the literature as the dimer 
(primer) model in case of a pair (triplet) of coupled 
wells] has been prompted recently by the construction 
of devices where Bose-Einstein condensates (BEC) inter- 
act through the tunneling effect (see 0] and references 
therein). The theoretic work focused on such models, 
both in the atomic physics community and in other ar- 
eas of theoretical physics, has supplied a large amount of 
results disclosing a quite structured interwell dynamics. 

The two-well model (TWM) -used to represent two 
coupled BECs in a symmetric double- well potential- has 
been investigated within a picture based on the alge- 
bra su(2) in Refs. ||, where, after stemming the model 
from the many-body quantum theory of BECs, the initial 
state with the atomic population self-trapped in one well 
is shown to evolve in delocalized oscillations involving 
both the wells. The same model has been studied previ- 
ously in Ref. j| , both at the quantum level and from the 
point of view of the dynamical system theory, to illustrate 
the level splitting that characterizes the dimer spectrum 
as a manifestation of the orbit bifurcation in the dimer 
phase space. In Refs. |g| the dynamics of the asymmetric 
TWM have been faced within the mean-field formulation 
relatively to the 7r-phase oscillations as well as the self- 
trapping effect. The latter was considered as well in Ref. 
Q and therein interpreted as a symmetry breaking phe- 
nomenon. More recently, the TWM (and its S-well gen- 
eralization) has been related ||j7]] to the Bose-Hubbard 
model [|| and the two- well ground-states have been inter- 
preted as insulator/superconducting regimes. In particu- 
lar, reformulating the TWM in an effective single-boson 
realization -generalizable to any S-well system- has been 
shown to favour the use of the system symmetries as well 
as the recognition of the inner parameters controlling the 
occurrence of doublets in the energy spectrum. 

In this paper we consider some recent results proving 
the existence of configurations with self-trapping within 
the dynamics of symmetric trimer (identical interwell 
couplings). These have been obtained in Ref. (I) by re- 



casting the trimer Hamiltonian within a two-boson op- 
erators picture (introduced in the sequel) which involves 
the algebra su(3). Such a picture is the extension of the 
dimer case || based on the su(2) (the formal setup for 
S-well models involves [||J^] the algebra su(S)). The main 
contribution of this paper is to apply to the trimer an 
alternative approach that both reproduces the results of 
the su(3) picture and show how the dynamical mecha- 
nism causing self-trapping not only depends on the tun- 
neling amplitude but also from the system initial con- 
ditions. Such an approach relies on a boson coherent 
state formulation previously developed for both boson 
and spin lattice models which seems to be very sim- 
ple and effective. The symmetric trimer is described by 
Hamiltonian 



H = UY J 



vN - 



(afa e + a\ai) 



with Tn — T, that one can derive from the many-body 
quantum theory of BECs through a three-mode expan- 
sion of the condensate field operator Q. Parameters U, 
v, T, account for the interatomic scattering, the exter- 
nal potential and the tunneling amplitude, respectively; 
rii = afai count the bosons in the ith well (N = Ejni), 
while the destruction (creation) operators a, (af) obey 
the canonical commutators [<2j,a^] = Sa. Preceding 
studies of the trimer dynamics have been focused on the 
asymmetric case characterized by tunneling amplitudes 
T12 > Ti3,T 2 3- Classically (aidf = afa,i, d\ = a|), the 
asymmetric trimer has revealed Jl0| the presence of ho- 
moclinic chaos, while, at the quantum level, the survival 
of breather configurations has been investigated on 
the trimer viewed as the smallest possible closed chain. 

If one derives the Heisenberg equations related to H for 
the boson operators Oj, af and implements the random 
phase approximation in the equations for their expecta- 
tion values Zi = (c^), z* — (af), the resulting equations 
for the three- well dynamics are (j = 1,2,3) 



ihij — (2U\zj 



+ T/2)z ] -T(z 1 +z 2 + z 3 )/2, (1) 



1 



which entail E conserved quantity replacing the 

total boson number N such that [TV, H] = 0. The Hamil- 
tonian structure of the Heisenberg equations is inherited 
by Eqs. (|l|) that, in fact, are also obtained from Ji(Z, Z*) 



T{z*z 3+ i + c.c.)/2], by using 



the standard Poisson brackets {z^, zj} — iS k j/h. 

Another significant way to obtain Eqs. ([j]) from H re- 
lies on applying the time-dependent variational principle 
on a suitable trial state | <J/) = e | Z) with Z = (z\ , z% , .. .) , 
where z r 's are time-dependent complex parameters ac- 
counting for the system evolution. Performing the varia- 
tion of (iff\(idt — H)\^>) = furnishes a system of hamil- 
tonian equations for Z = (zi, z r ) and identifies 9 with 
the action of the system. If the trial state is defined as M 



e i6 \zi) 



Z2 



*3 



(2) 



where \zi) are the standard bosonic coherent states that 
obey the defining equation Oi\zi) = Zi\zi), then Eqs. ([!]) 
are recoverd (up to the shift v — > v + U ) in which 
Zi = (zjla^Zi), z* = (zjla+lzi) = z*, \zi\ 2 = (zj\rii\zi), 
and d6/dt is the Lagrangian associated to H. In addition 
to describing the system evolution through | ty) , this ap- 
proach also provides a natural way to find the quantum 
configuration (in terms of states) corresponding to the 
initial conditions of a given classical motion. 

In Ref. jj| the semiclassical treatment of the trimer 
dynamics was based on deriving the equations of mo- 
tion for the expectation values of the two-boson op- 
erators forming the basis of su(3) instead of at, af . 
Such an algebra is generated by the creation operators 
t\ = ala 2 , t 2 = o-t a 3> e 3 = a 3~ a ii the destruction oper- 
ators ef — (ei) + , i — 1,2,3 and the (so-called) Cartan 
operators h 2 — [D 2 — D 3 )/^/3, hi = D\, where 



Di 



m — n 2 



Do 



n 2 - n 3 



Do 



n 3 - m 



(3) 
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By using imbalance operators (||), the su(3) algebraic 
structure is specified by the commutators [e^e/"] = 2Di, 
[ei,et] = Sitk^i, [De,e e ] = e e with i,k,£ € [1,3] (euu 
is the standard antisymmetric symbol), together with 
[ei,Dg\ — Ci/2, and [ei,e/] = 0, for i ^ I. Expressing 
Hamiltonian H through h\ and h 2 one finds 

2U{h\ + h 2 2 ) - f(N) - J(ei + e 2 + e 3 + h.c.) , (4) 



H 



with f(N) := UN 2 /3 + vN, where the operator TV is a 
group invariant, namely [N,g] = 0, Vg G su(3). This im- 
plies that [N, H] = 0. In this framework the Heisenberg 
equations are easily carried out. If the random phase ap- 
proximation (AB) = (A) (B) is also implemented Heisen- 
berg's equations for the su(3) generators take the form 



ie k = — (T + AUe k ) D k + ^e ku {4 ~ 4) . 



ih\ 
ih 2 



[(2 £l - eg - e 2 ) 
[\/3(e 3 - e 2 ) - 



c.c. 



(5) 



where we have used the displacement operators Dj for 
simplifying the formulas. Notice that, in Eqs. (|5|) the 
approximation (AB + BA) = 2(AB) has been repeatedly 
applied to bilinear terms, and q, ef, hi, h 2 have been 
used in place of their expectation values (q), (e^~), (hi), 
(h 2 ). A possible integrable regime is achieved by setting 

hi = (^ ni = n 2 ) , e 2 -e^ = 0, ei - = , 

which leads to the reduced system of equations 

= i(4 - £ 2) 

ie 2 = |(e 2 -e+-2 J D 2 )-4C/ J D 2 e 2 (6) 



ihi 



V3e 



c.c. 



Their solutions have been calculated implicitly by geo- 
metric arguments and reproduced numerically for various 
choice of initial conditions in Ref. [Q . 

In the alternative solution scheme based on Eqs. (jl]) 
the above constraints reduce to impose the condition 
Z\ = z 2 . This selects an integrable sub-dynamics. In 
fact, Eqs. ([T]) become two, 



ihzi = (2U\zi\ 2 - v)zi - |(z x + z 3 ) 
ihz 3 = (2U\z 3 \ 2 ~v)z 3 -Tzi , 



(7) 



where the two costants of motion corresponding to the 
energy and the total boson number (we set m = |zi| 2 ) 



E = U(2n\ + nf) - vN 
N = 2ni + n 3 



Tni - T(z 3 Z\ -I- z\z 3 ) 



(8) 



make Eqs. (0) integrable. The dynamical behavior is 
obtained explicitly via a standard quadrature procedure 
(see Refs. f|,|l2|) which furnishes the phase-independent 
equation for D 3 



Di 



— (4T 2 nm 3 
16 



R z 



(9) 



by substituting R := [E + vN + Tm - U(2n\ + n§)] 
= — T(zgZi + C.c.) inside the (squared) equation D 2 = 
— 9T 2 [zgZi — c.c. ] 2 /16 for D 3 . Introducing the further 
constant of motion N to obtain D 2 written in terms of 
the unique variable D 3 requires that ni and n 2 are ex- 
pressed as m = (N - 2D 3 )/3 and n 3 = [N + 4D 3 )/3. 
These, in turn, substituted in Eq. (H) give the equation 

D 2 = ^-(N-2D 3 )(N + 4D 3 ) -^-R 2 (D 3 ), (10) 
for the imbalance variable D 3 = (n 3 — ni)/2, in which 



T 



U 



c.c. 



R(D 3 ) =E + vN+ -(N - 2D 3 ) - -(N z + 8D 2 ) 



2 



= T) {(A 



D 2 )[T + W(A + D 2 )] - TNK(P)} 



(a + 2) 2 
9 3 (0)-6 



— 9a 2 ]2cosA, 
i(0). Thesec- 



with A := D 3 (0), K{P) := §[(. 
P := (a, A), a = 2A/JV, and A := 
ond version of R(D 3 ) is obtained by writing E in terms of 
the initial conditions D 3 (0), 6k(0). Phases 8j are defined 
by Zfc = Y / nfee l9fc . Eq. (|l0|) can be cast in the dimension- 
less form (dx/ds) 2 = -2V T (x; P) with s := NUt and 

V T {x-P) :=^[{a-x){a+x+T/2)-TK] 2 - T —{l-x){l+2x) 

where x := 2D 3 /N (x E [-1,1]), r := T/iVf/. In view 
of the fact that both the squared term in V T (namely 
R 2 ) and (dx/ds) 2 are nonnegative, the further condition 
(1 — x)(l + 2x) > must be accounted for which implies 
the restriction of the x range to — 1/2 < x < 1. 

The reduction of Eq. (0) to Eq. (|l^) allows one to 
construct explicit solutions in terms of elliptic func- 
tions by recasting the quartic term via standard trans- 
formation methods Jl4| . This will be enacted else- 
where. Operationally, our goal -the description of bifor- 
cation mechanism inherent in Eq. (]l0|)- can be achieved 
as well through the equivalent potential problem £ = 
^(dx/ds) 2 + V T (x; P) at £ = 0, where parameters N, K 
and x(0) in V T are fixed by setting the initial conditions. 

With negative r and a suitable choice of the other pa- 
rameters, V T can exhibit an asymmetric double- well. In 
general, three solutions are obtained by annihilating 



dV T „ o 3r 9 1 

— - = 2x 3 H x 2 + - 

dx 2 4 



9r 2 -8/3 T (P) 



T 

X ~2 



where /3(r, P) := (a + t/A) 2 - t(K + r/16), that corre- 
spond to a maximum of V T (x; P) with two side minima. 

In particular, setting a = 1 reproduces the conditions 
under which dynamics was studied in Ref. (depleted 
twin wells, that is ^3(0) = 1), and leads to the potential 

1 2 2 

V T (x) = -(1 - x) 2 (x + T - + l) - T -(\ - x){\ + 2x) 

whose maximum is such that V T (x m ) = with x m = 
when t = —2/3. For r > —2/3 one has V T (x m ) > 0. 
The important feature thus emerging (see Fig. 1) is that, 
whenever the potential maximum is nonnegative, V T (x) 
generates two noncommunicating basins with V T (x) < 
(separated by a forbidden interval where V T > 0) entail- 
ing two independent oscillatory motions. In each basins 
the motion has a periodic character. This represents the 
bifurcation effect reminescent of the behavior manifested 
by two- well dynamics |I],[f| . 

What we emphasize here, based on our Zj description, 
is that the onset of separated motions can be caused by 
varying the other parameters of the problem. In particu- 
lar, a high sensitivity is manifested relatively to the initial 
phases incorporated in A. Suitable changes of the latter 



are capable of switching on the bifurcation mechanism 
even for o / 1. Such a situation is represented in Fig. 2 
for a — 0.99 (twin wells almost empty) and r = —2/3, 
where various potential wells are generated by varying 
cosA in [—1,1]. For sufficiently low values of cosA the 
presence of the maximum is ensured. The 'opposite' case 
a = —0.49 and r = —1/3 (corresponding to twin wells 
almost half- filled and n 3 (0) ~ 0) of Fig. 3 confirms the 
presence of isolated basins as well as the case with a more 
negative coupling r = —0.7 < —2/3 and a — 0.99. 

Decreasing sufficiently the value of r (Fig. 4 illustrates 
the case r = —0.8 with a = 0.99) by keeping the same 
range for cosA entails situations where wells never ex- 
hibit a local maximum. This can be proved analitically 
in the special case r = — 1 in which the potential becomes 



Vr(x;P) = l 



(a-1/4) 2 + K - X< 



- + X 2 
16 



with X = x— 1/4, and the stationary points can be calcu- 
lated explicitly. One finds a maximum at x m = 1/4 with 
V T (x m ) < so that no bifurcation effect occurs. The side 
minima are placed at x T ^i — l/4± [A'— 1 + (a— 1/4) 2 ] 1 / 2 . 
These are real provided K — 1 + (a — 1/4) 2 > namely if 

cosA > [1 - (a - l/4) 2 ]/[(l - a)(l + 2a)]* . 

For a generic r, the maximum depends on a and A 
in a complicated way which makes difficult the analytic 
calculation of V T (x m ) and of its sign. Nevertheless, some 
necessary conditions ensuring its existence can be ob- 
tained explicitly. As suggested by Figs. 1-3, increasing A 
with both t and a constant implies that the maximum at 
x = x rn and the left minimum at x — xi reach the (flex) 
point x = c for critical value A = A*. Since the interval 
[xe, x m ] where dV T /dx > vanishes for xe, x rn — > c then 

tim m (dV r /dx) xijXm = = (d 2 V T /dx 2 ) c . (11) 

The derivation of the roots of d 2 V T /dx 2 — at x = c 

x± = ^{-1 ± [8(2a 2 + ar - 2Ar)/(3r 2 ) - 5] 1/2 } , 

from d 2 V T /dx 2 = 6x 2 + 3rx- 2/3 T (P) + 9r 2 /4 allows one 
to exploit the fact that the lowest one, x_, is a maximum 
of dV T /dx corresponding to the V T flex point at x = c. 
When 



(dV T /da 



3 3 / 



Nr [/3 T (P)]l-r(r+l)>0 (12) 



becomes negative the maximum disappears (see, e. g., 
Figs. 1-3). The bifurcation condition V T (x) > must be 
searched within the parameter space domain where a, A, 
t satisfy formula (p^). 

The analysis just developed shows that changing A 
can modify deeply the system dynamics and that, in gen- 
eral, the onset of the bifurcation effect is governed by the 



3 



complex interplay of all parameters a, A, r. A complete 
study of the dynamics requires that one considers any 
possible initial condition for the dynamics and thus, e. 
g., the situations in which ni(0) ^ ^2(0), excluded in the 
present paper. In this case the nonintegrable character 
of the system crops up in a dramatic way as illustrated 
in Fig. 5. The systematic analysis of fixed points for the 
symmetric three-well dynamics and thus the emergence 
of chaos close to the hyperbolic points is in progress at 
this moment. It will be discussed in a separate paper. 
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FIG. 1. By varying r in [—0.75, —0.63] with a = 1, V T (x,P) generates a second (small) basin on the left (dashed potential 
corresponds to r ~ —0.66). 



FIG. 2. By varying A in [0, ir] with a — 0.99, V T (x, P) generates a second (small) basin on the left (dashed potential 
corresponds to A ~ 1.40). 

FIG. 3. Representation of bifurcation mechanism by varying A £ [0, n] in V T (x, P) with a = —0.49, r = —1/3. 



FIG. 4. V T (x, P) with a — 0.99 and r = 0.8: a sufficiently negative r involves a single basin for any A £ [0, w]. 

FIG. 5. Poincare section of the £i — (j>\ plane, where £i := 1 — 2ni/N and 4>i =62 — 61, obtained by setting n 3 ~ 6.85; this 
is derived by numerical integration of Eqs. (1) with energy E ~ 92.33, N = 10, T = U = 1. 
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